Intro

Fit Tweedie models to biomass density of cod, flounder, plaice and dab (juveniles and adults) between 1993-2020 in the Baltic Sea using sdmMTB fit different oxygen and temperature-related covariates. Predict on grid and use those predictions to calculate trends and velocities.

Load packages & source functions

# Load libraries, install if needed
library(tidyverse)
#> Warning: package 'tidyr' was built under R version 4.0.5
library(readxl)
library(tidylog)
library(RCurl)
library(devtools)
library(sdmTMB)
library(rMR)
library(kableExtra)

# Source code for plots
source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
#> Warning: package 'sf' was built under R version 4.0.5

theme_set(theme_plot())
# Read data
d <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>%
  rename(X = x, Y = y) %>% 
  drop_na(depth) %>% 
  mutate(depth_ct = depth - mean(depth),
         depth_sc = depth_ct / sd(depth),
         depth_sq = depth_sc^2,
         year_f = as.factor(year),
         quarter_f = as.factor(quarter)) %>% 
  pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
                 flounder_juvenile, plaice_adult, plaice_juvenile),
               names_to = "group", values_to = "density") %>% 
  separate(group, into = c("species", "life_stage"), remove = FALSE)
#> Rows: 9376 Columns: 23
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (3): country, haul_id, ices_rect
#> dbl (20): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed 2 variables (X, Y)
#> 
#> drop_na: removed 112 rows (1%), 9,264 rows remaining
#> 
#> mutate: new variable 'depth_ct' (double) with 3,707 unique values and 0% NA
#> 
#>         new variable 'depth_sc' (double) with 3,707 unique values and 0% NA
#> 
#>         new variable 'depth_sq' (double) with 3,707 unique values and 0% NA
#> 
#>         new variable 'year_f' (factor) with 28 unique values and 0% NA
#> 
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> 
#> pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult, …) into (group, density) [was 9264x28, now 74112x22]

# Read metabolic parameter estimates and left_join
# mi_pars <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/mi_params.csv") %>% 
#   dplyr::select(-...1, -temp)

# TODO: for now we'll use plaice parameters for flounder, see "00_estimate_mi_params.Rmd"
# mi_pars2 <- mi_pars
# mi_pars2[4, c(4:5)] <- mi_pars2[1, c(4:5)]
# 
# d <- left_join(d, mi_pars2, by = "species")

# Read size csv to calculate the metabolic index
# sizes <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/sizes.csv") %>% 
#   mutate(group = paste(species, name, sep = "_")) %>% 
#   dplyr::select(group, B)
# 
# d <- left_join(d, sizes, by = "group")

Read the prediction grid

# Read prediction grid
pred_grid <- bind_rows(readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/master/data/clean/pred_grid_(1_2).csv"),
                       readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/master/data/clean/pred_grid_(2_2).csv"))
#> Rows: 454412 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): ices_rect
#> dbl (11): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Rows: 454412 Columns: 12
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (1): ices_rect
#> dbl (11): X, Y, year, lon, lat, depth, quarter, oxy, temp, sal, sub_div
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

Calculate the metabolic index

# # Before we can calculate the pressure based metabolic index, we'll need to convert oxygen from ml/L mg/L and partial pressure
# # Some conversions
# # oxygen: https://www.ices.dk/data/tools/Pages/Unit-conversions.aspx
# # pressure: https://bluerobotics.com/learn/pressure-depth-calculator/
# # First calculate pressure in pascal: p_tot = (r*g*h) + p_atm, then to kgPa and finally atm
# rho = 1030 # density of water
# p_atm = 1.01325 # atmospheric pressure
# g = 9.807 # gravity of earth
# 
# # Convert to partial pressure oxygen from ml/L
# d$po2 <- DO.unit.convert(d$oxy*(1/0.7), # d$oxy is in unit ml/L so we first need to convert to mg/L
#                          DO.units.in = "mg/L", 
#                          DO.units.out = "PP", 
#                          bar.press = ((rho*g*d$depth + p_atm)*10^-5)*0.986923, 
#                          bar.units.in = "atm",
#                          bar.units.out = "kpa",
#                          temp.C = d$temp,
#                          salinity = 10, # TODO: d$sal, # let's fix salinity for easier comparison. It has a minor effect in this range
#                          salinity.units = "pp.thou")
# 
# # Oxygen is ml/L, We want micro mol/L. 1 ml/l = 10^3/22.391 = 44.661 micro mol/l
# d$oxy_si <- (d$oxy * (10^3)) / 22.391
# 
# # Calculate metabolic indices for a given mass and the temperature and oxygen in data
# # Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
# kb <- 0.000086173324 # Boltzmann's constant
# t_ref <- 15 # arbitrary reference temperature
# 
# # Calculate the metabolic index for pressure and concentration 
# d <- d %>%
#   mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
#          phi = A0_po2*(B^n_po2)*po2 * exp((E_po2/kb)*inv_temp),
#          phi_c = A0_o2*(B^n_o2)*oxy_si * exp((E_o2/kb)*inv_temp))

Fit models and save AIC

# Fit models and save grid-level predictions
data_list <- list()

for(i in unique(d$group)) { 
    
    dd <- d %>%
      filter(group == i) %>%
      droplevels() %>% 
      drop_na(density)
    
    mesh <- make_mesh(dd, c("X", "Y"), n_knots = 200)

    plot(mesh)
    
    m <- sdmTMB(
      density ~ 0 + year_f + quarter_f + depth_sc + depth_sq,
      data = dd, mesh = mesh, family = tweedie(link = "log"),
      spatial = "on", time = "year", spatiotemporal = "IID",
      control = sdmTMBcontrol(newton_loops = 1))

    sanity(m)

    # Plot residuals
    res <- residuals(m) %>% as.data.frame()
    ggplot(res, aes(sample = V1)) + stat_qq() + stat_qq_line() + theme(aspect.ratio = 1)
    ggsave(paste0("figures/supp/qq_density_sdm/sdm_02_", i, ".pdf"),
           width = 17, height = 17, units = "cm")
    
    # Save model object
    saveRDS(m, paste("output/models/sdm_02_", i, ".rds", sep = ""))
    
    # Predict on grid
    pred_grid <- pred_grid %>%
      mutate(year_f = as.factor(year),
             depth_ct = depth - mean(dd$depth),
             depth_sc = depth_ct / sd(dd$depth),
             depth_sq = depth_sc^2,
             quarter_f = as.factor(quarter))

    pred <- predict(m, newdata = pred_grid) %>% mutate(model = i)
    
    data_list[[i]] <- pred
    
}
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 13 rows (<1%), 9,251 rows remaining
#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: new variable 'year_f' (factor) with 28 unique values and 0% NA
#>         new variable 'depth_ct' (double) with 11,144 unique values and 0% NA
#>         new variable 'depth_sc' (double) with 11,144 unique values and 0% NA
#>         new variable 'depth_sq' (double) with 11,144 unique values and 0% NA
#>         new variable 'quarter_f' (factor) with 2 unique values and 0% NA
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 13 rows (<1%), 9,251 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 629 rows (7%), 8,635 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 629 rows (7%), 8,635 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✖ `b_j` standard error may be large
#> ℹ Try simplifying the model, adjusting the mesh, or adding priors
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 235 rows (3%), 9,029 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 235 rows (3%), 9,029 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 622 rows (7%), 8,642 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA
#> filter: removed 64,848 rows (88%), 9,264 rows remaining
#> drop_na: removed 622 rows (7%), 8,642 rows remaining

#> ✔ Non-linear minimizer suggests successful convergence
#> ✔ Hessian matrix is positive definite
#> ✔ No extreme or very small eigenvalues detected
#> ✔ No gradients with respect to fixed effects are >= 0.001
#> ✔ No fixed-effect standard errors are NA
#> ✔ No standard errors look unreasonably large
#> ✔ No sigma parameters are < 0.01
#> ✔ No sigma parameters are > 100
#> ✔ Range parameter doesn't look unreasonably large
#> mutate: changed 908,824 values (100%) of 'depth_ct' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sc' (0 new NA)
#>         changed 908,824 values (100%) of 'depth_sq' (0 new NA)
#> mutate: new variable 'model' (character) with one unique value and 0% NA


# Save predictions and sims as data frames
df <- dplyr::bind_rows(data_list)

write_csv(df, "output/grid_pred_biomass_02.csv")